Vibrational Corrections to NMR Spin–Spin Coupling Constants from Relativistic Four-Component DFT Calculations

Zero-point vibrational (ZPV) corrections to the nuclear spin–spin coupling constants have been calculated using four-component Dirac–Kohn–Sham DFT for H2X (where X = O, S, Se, Te, Po), XH3 (where X = N, P, As, Sb, Bi), and XH4 (where X = C, Si, Ge, Sn, and Pb) molecules and for HC≡CPbH3. The main goal was to study the influence of relativistic effects on the ZPV corrections and thus results calculated at relativistic and nonrelativistic approaches have been compared. The effects of relativity become notable for the ZPV corrections to the spin–spin coupling constants for compounds with lighter elements (selenium and germanium) than for the spin–spin coupling constants themselves. In the case of molecules containing heavier atoms, for instance BiH3 and PbH4, relativistic effects play a crucial role on the results and approximating ZPV corrections by the nonrelativistic results may lead to larger errors than omitting ZPV corrections altogether.


■ INTRODUCTION
The standard approach to calculations of molecular properties within the Born−Oppenheimer approximation is to evaluate them at some reference geometry, usually the equilibrium geometry. However, it is well known that high-precision calculations of molecular properties require taking into account vibrational corrections. 1,2 This is particularly true of the NMR properties: nuclear spin−spin coupling constants and nuclear shielding constants, which both are sensitive to geometry distortions and thus to effects associated with nuclear motion.
There are several approaches for evaluating vibrational corrections to the spin−spin coupling constants, 3−5 differing in accuracy and computational cost. The majority of the effect can be approximated by computing the zero-point vibrational (ZPV) corrections, 5 that is, the difference between the equilibrium value and the averaged value for the ground vibrational state. ZPV corrections are usually calculated by perturbation theory 6−8 and included in accurate computational studies.
On the other hand, it is well known that relativistic effects (understood as a difference between the results obtained using relativistic and nonrelativistic Hamiltonians) on NMR parameters can be non-negligible already for third-row elements. 9 When both relativistic and vibrational corrections need to be accounted for, it is usually done by an incremental approach: calculating zero-point vibrational corrections using a nonrelativistic Hamiltonian and adding them to the relativistic value. This assumes that the nonrelativistic property and energy surfaces are sufficiently close to being parallel to the correct relativistic ones, or, in other words, that the relativistic corrections are similar for all geometries close to the equilibrium geometry. For many systems, this approach has been applied successfully 10,11 but it is not always the case: it has been shown that in some cases 12 derivatives of the spin− spin coupling constants with respect to internuclear distance can even differ in sign when calculated with nonrelativistic and relativistic Hamiltonians. There is, therefore, a need to calculate also ZPV corrections at the relativistic level of theory in order ensure correct estimates for these effects.
■ METHODS Theory. The most popular approach to calculating vibrational corrections to NMR parameters is the approach of Kern et al., 6−8 in which second-order perturbation theory is used. It has also been applied in the present work. It should be noted that this method implies only small-amplitude nuclear motions. In the case of large-amplitude nuclear motions (e.g., internal rotation) other methods, for example, molecular dynamics, must be employed, 13−15 as it is important to distinguish conformational equilibria from large-amplitude motions.
In the perturbational approach, the unperturbed groundstate vibrational wavefunction is written as a product of harmonic oscillator wavefunctions in normal coordinates 5,16 where ϕ K n (Q K ) is the nth excited harmonic oscillator state of the Kth normal vibrational mode, and the summation runs over 3N − 6 normal modes, where N is the number of atoms in the molecule. In the next step, a full set of virtual excitations from Ψ (0) (Q) is used to expand the first-order correction to the ground-state vibrational wavefunction, Ψ (1) (Q). If the formula for Ψ (1) (Q) is limited to the third-order Taylor expansion of the potential energy surface, the only relevant contributions are from single and triple excitations Φ K 1 (Q) and Here, for example, Φ KLM ABC (Q) has been obtained from Φ 0 (Q) by exciting the Kth, Lth, and Mth modes to the Ath, Bth, and Cth harmonic oscillator states, respectively. The expansion coefficients in the above can be written (in atomic units) as 17 i where and ω K is the mass-weighted harmonic frequency for the Kth normal mode. In the equilibrium geometry F K = 0. A vibrationally averaged molecular property P can be now calculated as an expectation value If P is expanded in a Taylor series about the equilibrium geometry The final form of the formula for the ZPV correction to a property P is therefore The first term in the above equation is the harmonic contribution to the ZPV correction and the second term is the anharmonic contribution.
This formula has been used in the present work to calculate ZPV corrections to the nuclear spin−spin coupling constants.
Implementation. Our program works as an external driver to the Dirac 18 program package, but can in principle be adapted to any other program. The ZPV corrections to the spin−spin coupling constants are calculated with the approach of Kern et al. 6−8 using eq 13. In the case of NMR parameters, there is no analytic implementation for the energy and property derivatives and thus the method is fully numerical, which means that the first and diagonal second derivatives of the spin−spin coupling constants, as well as the harmonic frequencies and the semi-diagonal part of the cubic force field, are calculated numerically.
Numerical Derivatives. The molecular Hessian, normal coordinates, and vibrational frequencies are calculated as described in our previous paper. 19 Once the vibrational frequencies and normal coordinates are computed, the first and second derivatives of the spin−spin coupling constants with respect to geometric distortions along the normal The approach thus involves performing a number of energy and property calculations, in which atoms are being displaced from their original positions along the normal coordinates. In the case of a nonlinear N-atom molecule (with 3N-6 vibrational modes), 45N 2 − 165N + 150 single-energy computations and 6N − 11 property computations need to be run to determine the ZPV corrections.
When carrying out numerical differentiation, it is essential that an appropriate step length (h in the above equation) is used to ensure numerically accurate results. On one hand, if the step length is too small, numerical errors will dominate due to the approximate solution of the perturbed wavefunctions. On the other hand, if it is too large, the derivatives will be contaminated by higher-order terms. We have performed test calculations of the ZPV correction for the water molecule with a number of different step lengths in the range of 0.1 bohr amu 100 bohr amu · · . The calculations turned out to be numerically stable for step lengths between 1 bohr amu and 50 bohr amu · · . Based on the above, for all subsequent calculations we have used a step length of 10 bohr amu · .
■ COMPUTATIONAL DETAILS Geometry Optimization. Geometry optimizations have been performed using the Dirac 18 program at the same level of theory as the ZPV correction calculations carried out afterward in order to ensure that the molecular gradient is zero (a condition for the harmonic approximation). The convergence threshold for the gradient was 10 −4 au. For comparison, also nonrelativistic calculations have been carried out. In the case of the nonrelativistic computations, the speed of light has been scaled to 2000.0 au in the Dirac− Coulomb Hamiltonian.

Single-Point Energy and Property
Because the semi-diagonal part of the cubic force field was calculated numerically, the convergence threshold for all the single-point energy calculations needed to be tight. For this reason, the convergence threshold for the error vector was set to be 10 −10 and in a few cases (about 10%) 10 −8 if the number of iterations exceeded 50.
Molecules under Investigation. In order to test the newly developed method for calculating ZPV corrections to spin−spin coupling constants, simple systems consisting of 3, 4, and 5 atoms have been chosen: • H 2 X where X = O, S, Se, Te, Po; • XH 3 where X = N, P, As, Sb, Bi; and • XH 4 where X = C, Si, Sn, and Pb.
For some of these systems, vibrational corrections to the nuclear spin−spin coupling constants are known in the literature. 29−31 In addition to this, to illustrate the usefulness of the method for larger systems, we have calculated ZPV corrections to the spin−spin coupling constants for an acetylene derivative, HC�CPbH 3 .
As the vibrational frequencies are incorporated in the formula for the ZPV correction (see eq 13) and vibrational frequencies change for different isotopes of the same element, we needed to select the isotopic constitution of the molecules for which the calculations were performed. In the case of J(H− X) couplings, 1 H and the most abundant magnetic isotopes of element X ( 17 O, 33 S, 77 Se, 125 Te, 209 Po, 14 N, 31 P, 35 As, 123 Sb, 209 Bi, 13 C, 29 Si, 73 Ge, 119 Sn, and 207 Pb) were chosen (although we are aware that for many of them, the measurements of the spin−spin coupling constants are not possible because of the quadrupole moment of the nucleus and thus the associated line broadening). In the case of J(H−H) couplings, the computations were carried out for 1 H and the most abundant isotope of element X: 16 O, 32 S, 80 Se, 130 Te, 209 Po, 14 N, 31 P, 35 As, 121 Sb, 209 Bi, 12 C, 28 Si, 74 Ge, 120 Sn, and 207 Pb. As far as the HC �CPbH 3 molecule is concerned, in order to limit the computational cost, the calculations were run only for 1 H, 13 C, and 207 Pb.
■ RESULTS AND DISCUSSION Spin−Spin Coupling Constants. Even though the main focus of this work is to analyze the role that relativistic effects play on the ZPV corrections to spin−spin coupling constants, the results for the spin−spin coupling constants themselves will be briefly discussed for the sake of completeness. They have been collected in Table 1.
In the case of couplings that involve the X atoms, which have different magnetogyric constants, we discuss reduced spin− spin coupling constants, K, due to their independence with the magnetogyric constants. Relativistic effects are noticeable and All of the above findings are in line with previous studies. 35 −38 Effects of Relativity on the First and Second Derivatives of Spin−Spin Coupling Constants. The ZPV corrections to the spin−spin coupling constants depend on the first and second derivatives of the coupling constants with respect to nuclear distortions, the cubic force field, and the harmonic vibrational frequencies. Each of these parameters can to a different extent be sensitive to relativistic effects. We have, therefore, also investigated the influence of relativity on the first and second derivatives of the coupling constants with respect to normal coordinates. The results calculated with the  relativistic and nonrelativistic approaches are shown in Table 2, for the sake of brevity only for the H 2 X systems. All the following observations can be generalized to the XH 3 and XH 4 systems. Analysis of the relativistic and nonrelativistic results in Table  2 indicates that the relativistic effects tend to be more pronounced for the derivatives of the coupling constants than for the coupling constants themselves. This is true both for the first and second derivatives of the coupling constants.
When analyzing the results, two interesting observations can be made. First of all, the derivatives with respect to different normal coordinates show different sensitivity to relativity. For instance, in the case of H 2 Te, the relativistic value for ).

ZPV Corrections to Spin−Spin Coupling Constants.
The results of the calculations of ZPV corrections to the spin− spin coupling constants computed with both relativistic and nonrelativistic methods for H 2 X, XH 3 , and XH 4 are presented in Tables 3 and 4 for 1 K XH and 2 J HH , respectively.
Because a method for calculating ZPV corrections to NMR parameters is implemented in the Dalton 39,40 program, some nonrelativistic calculations have been performed with this program in order to check the consistency of the approach. All of the Dalton computations have been run with the same uncontracted basis set and exchange−correlation functional as above. The results can be found in the Supporting Information. In almost all cases, Dalton produces results that are in excellent agreement with the results obtained with our newly implemented method. The only exception is the ZPV correction to 1 J TeH , for which the result obtained with Dalton is unphysically large, suggesting a problem with this calculation.
Effects of Relativity on 1 K XH ZPV . As shown in Table 3, relativistic effects to the ZPV corrections of 1 K XH become noticeable for lighter systems than was the case for the spin− spin coupling constants themselves. For H 2 Se, PH 3 , AsH 3 , and GeH 4 , the differences between nonrelativistic and relativistic results for the total ZPV correction fall within the range of 10− 15% of the relativistic value.
The most striking differences between the ZPV corrections to 1 K XH calculated with nonrelativistic and relativistic approaches occur for SbH 3 , BiH 3 , and PbH 4 . In the case of SbH 3 and BiH 3 , 1 K XH ZPV changes from −14.01 to −2.52 × 10 19 · m −2 ·kg·s −2 ·Å −2 and from −18.35 to −143.30 × 10 19 ·m −2 ·kg· s −2 ·Å −2 , respectively. We note that an observed decrease or increase in the value is the same for the spin−spin coupling constant and the corresponding ZPV correction when the method is changed from nonrelativistic to relativistic. The nonrelativistic absolute value of the ZPV correction to the coupling constant for SbH 3 is larger than the relativistic value of the coupling constant itself, whereas the relativistic value of the ZPV correction constitutes about 20% of the relativistic value of the coupling constant.
An interesting observation can be made for PbH 4 . As the spin−spin coupling constants increase significantly using a relativistic Hamiltonian, the ZPV correction decreases by almost 150%. Furthermore, the nonrelativistic ZPV correction constitutes around 33% of the nonrelativistic coupling constants, whereas this percentage decreases to only 5% for the relativistic results.
In almost all cases, a change in the method from nonrelativistic to relativistic leads to changes in both the harmonic and anharmonic terms that are mostly of the same magnitude, with the two notable exceptions of H 2 Po and PbH 4 . For H 2 Po, the change in the harmonic term is 127%, whereas the change in the anharmonic term is 34%, and for PbH 4 these changes are 488 and 91%, respectively.
As our main goal is to study relativistic effects on ZPV corrections to spin−spin coupling constants rather than reproduce experimental results, experimental values were not given in Tables 3 and 4. A brief comparison with experimental data in gas phase 41 and vibrationally averaged reduced spin− spin coupling constants, 1 K XH , calculated at the relativistic level is given in Table 5 for CH 4 , SiH 4 , GeH 4 , and SnH 4 . It is clear that in the case of CH 4 and SiH 4 , adding the ZPV correction does not bring the calculated spin−spin coupling constants  As for 2 J HH , in almost all cases the relative change in the harmonic and anharmonic terms is of the same magnitude when nonrelativistic and relativistic results are compared, the only exceptions being H 2 Po, AsH 3 , and SnH 4 .
Effects of Relativity on ZPV Corrections to Spin−Spin Coupling Constants for HC�CPbH 3 . The results of calculations of spin−spin coupling constants and the corresponding ZPV corrections for HC�PbH 3 computed with both relativistic and nonrelativistic methods are given in Table 6. The results are also compared to experimental values.
As far as the comparison of relativistic and nonrelativistic values of the spin−spin coupling constants is concerned, not surprisingly, relativistic effects play a key role in the case of 1 J CPb , 2 J CPb , and 3 J PbH . The HALA effect is almost non-existent for 1 J HC and 1 J CC , whereas 2 J HC (geminal coupling with the Pb atom in the middle) decreases by over 10% when a relativistic approach is used.
Using a relativistic Hamiltonian in the calculations of ZPV corrections turns out to be important both for spin−spin coupling constants that involve and do not involve a heavy atom. Relativistic effects constitute from 6% (for 1 J HC ) to as much as 297% (for 3 J PbH ) of the total relativistic ZPV correction. An interesting observation can be made for the ZPV correction to 2 J CPb . Even though the differences between the total ZPV corrections calculated with relativistic and nonrelativistic methods are relatively small, the changes of harmonic and anharmonic contributions are much larger. The harmonic contribution increases and the anharmonic con-tribution decreases and these changes partially cancel each other in the total value of the ZPV correction. The cancellation of the relativistic effect is thus coincidental, and in other cases, the ZPV corrections on the one-bond couplings of this type may be much more affected by relativity, as seen for the H 2 X, XH 3 and XH 4 systems.
The available experimental data refer to the ethylenesubstituted acetylene derivative HC�CPb(C 2 H 5 ) 3 , whereas the coupling constants and ZPV corrections discussed below have been calculated for compounds containing hydrogen atoms instead of ethylene groups. In ref 12, the influence of such a substitution was studied and a correction to the experimental value for HC�CPb(C 2 H 5 ) 3 can be introduced so as to estimate an "experimental" value for HC�CPbH 3 . These values are given in parentheses next to the experimental values for HC�CPb(C 2 H 5 ) 3 in Table 6. It can be noticed that for 1 J HC , 2 J HC and 1 J CC , adding the ZPV calculated using a relativistic approach brings the spin−spin coupling constants closer to the estimated "experimental" value, whereas for 1 J CPb and 2 J CPb the ZPV correction brings the calculated coupling constant further from the estimated "experimental" value. However, the vibrational effects are not the only effects that should be taken into account when comparing computational results to experiment. A study of available experimental data shows that in this case, solvent effects might also play an important role. 42 Moreover, the remaining disagreement with experiment might also be due to the errors resulting from the use of DFT with the B3LYP functional.

■ CONCLUSIONS
We have presented a numerical method for calculating the ZPV corrections to spin−spin coupling constants with relativistic four-component DFT. Test calculations have been performed for hydrides of elements from groups 14, 15, and 16, and for HC�CPbH 3 in order to demonstrate the versatility of the method.
For both the ZPV corrections to spin−spin coupling constants and the derivatives of the spin−spin coupling constants, the effects of relativity become notable much earlier in terms of the atomic number of the heavy element, for example selenium and germanium, compared to the spin−spin coupling constants. Moreover, our calculations demonstrate that as far as molecules containing heavier atoms are concerned, for instance BiH 3 and PbH 4 , relativistic effects have such a great impact on the results that the commonly In addition to this, ZPV corrections to spin−spin coupling constants have been computed for HC�CPbH 3 . Relativistic effects turned out to be at least noticeable, if not crucial, for all the calculated ZPV corrections to spin−spin coupling constants. Analysis of the results obtained shows that relativity should be taken into account for couplings that involve a heavy atom.
■ ASSOCIATED CONTENT
Comparison of results for ZPV corrections to spin−spin coupling constants calculated with Dalton and our newly implemented method and different optimized geometries (PDF)